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A quasi-diagonal approach to the estimation of Lyapunov spectra 
for spatio-temporal systems from multivariate time series 
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We describe methods of estimating the entire Lyapunov spectrum of a spatially extended system 
from multivariate time-series observations. Provided that the coupling in the system is short range, 
the Jacobian has a banded structure and can be estimated using spatially localised reconstructions 
in low embedding dimensions. This circumvents the "curse of dimensionality" that prevents the 
accurate reconstruction of high-dimensional dynamics from observed time series. The technique is 
illustrated using coupled map lattices as prototype models for spatio-temporal chaos and is found 
to work even when the coupling is not strictly local but only exponentially decaying. 
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I. INTRODUCTION 

One of the most important tools for investigating 
chaotic dynamical systems is the spectrum of Lyapunov 
exponents. These exponents measure the asymptotic ex- 
ponential divergence or convergence of two infinitesimally 
close orbits. In this paper we are interested in estimating 
all of the Lyapunov exponents of a spatio-temporal sys- 
tem from an observed multivariate time-series without 
prior knowledge of the dynamics governing the system. 
Hitherto, most efforts in this area have concentrated on 
only estimating the largest (or few largest) exponent(s). 
However, at least some of the negative exponents are 
needed if we wish to estimate the dimension of the at- 
tractor via the Kaplan- Yorke conjecture. Estimating all 



of the Lyapunov exponents from a time series for spa- 
tially extended systems is a daunting task. The funda- 
mental problem is the high dimensionality of the system 
which prevents an accurate reconstruction of the dynam- 
ics from observed data. In particular, whatever method 
of function approximation we use to reconstruct the dy- 
namics, we need to have a reasonable spread of data in 
the region of interest. Thus, for instance, a local linear 
or quadratic approach estimates the value of a function 
at a point y £ R d by performing a linear least squares re- 
gression in a neighbourhood of y. Such a neighbourhood 
must contain a sufficient number of data points to yield 
a meaningful estimate and yet not be so large that the 
function is no longer linear (or quadratic respectively). 
As d grows, more and more data is necessary to ensure 
that sufficient numbers of neighbours can be found and 
for d much larger than 6 or so the amount of data required 
becomes completely impracticable. This problem is often 
informally referred to as the "curse of dimensionality" . 

It is therefore perhaps surprising that, at least in some 
cases, the whole Lyapunov spectrum can be estimated 
successfully from observed data, see Ref. [jl| and || . Both 
of these papers focus on a lattice of locally coupled fully 
chaotic logistic maps, though Ref. J| (here referred as 
ORS99) also considers the effects of skewing the map (as 
we also do below in VI), and of replacing it by a much 
more non-linear function. The results for the latter two 
cases are substantially less satisfactory than for the stan- 
dard logistic map. It turns out that the key property 
that makes it possible to obtain reasonable estimates of 
the spectrum for the logistic map lattice using the meth- 
ods in ORS99 is that the dynamics at one spatial location 
is easily approximated within the space of functions used 
to fit the dynamics. Hence what at first sight appears 
to be a local quadratic fit is in fact a global fit. This al- 
lows the "curse of dimensionality" to be circumvented, 
and good estimates of the dynamics to be obtained even 
in high dimensions. As the local dynamics moves away 
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from the space of functions used to fit the dynamics the 
estimates of the Lyapunov spectrum rapidly deteriorate. 
One approach to overcoming this might be to attempt to 
estimate the local dynamics and then use a suitable basis 
to fit the dynamics Q. 

The alternative, which is the method used in Ref. 
(here referred as BH99) and which we pursue here, is 
based on the observation that as long as the spatial cou- 
pling in the system is reasonably short range, optimal 
predictions are often obtained in very low embedding di- 
mensions Q , even when the dimensionality of the attrac- 
tor is high. This somewhat counter-intuitive result may 
be explained by the rapid spatial decay of the dependence 
of the dynamics at one site on its neighbours ||. This 
suggests that it ought to be possible to obtain reasonable 
estimates of the Jacobian by performing appropriate fits 
in low embedding dimensions, hence avoiding the "curse 
of dimensionality" . 

The aim of this manuscript is to describe a method 
based on this intuitive idea and evaluate its performance 
in various circumstances. The method is based on only 
reconstructing the non-zero entries of the Jacobian in a 
row-by-row fashion. More precisely if the coupling in a 
spatio-temporal system is reasonably local then the only 
non-zero entries in the Jacobian are near the diagonal. 
Rather than estimating the whole Jacobian, one should 
only estimate the non-zero entries (which can be done in a 
low embedding dimension) and set the remaining entries 
to 0. We call such an approach a quasi-diagonal one. 
Essentially the same technique is used by Biinner and 
Hegger in BH99 who successfully apply it to the logistic 
coupled map lattice, using local linear fits. However, as 
ORS99 shows, the standard logistic lattice is a relatively 
easy system for which to estimate the spectrum, even 
using local linear fits (rather than quadratic ones, as in 
ORS99 and here), ft is thus impossible to judge from 
BH99 how well a quasi-diagonal approach works when 
the local dynamics presents a bigger challenge, or indeed 
when the coupling is anything but nearest neighbour. Fi- 
nally, ORS99 demonstrates that dramatically improved 
estimates of the spectrum can be obtained by truncat- 
ing the outer layer(s) of the estimated Jacobian, thereby 
eliminating boundary effects. Such truncation is not con- 
sidered in BH99. 

The present work combines the best features of BH99 
and ORS99 together with additional generalisations and 
extensive supporting numerical evidence. In particu- 
lar, we evaluate the quasi-diagonal approach when ap- 
plied to a more difficult local map, using a more flex- 
ible (quadratic) fitting basis. We investigate the effect 
of truncating the outer layers of the estimated Jacobian 
and assess the effectiveness of our approach when the cou- 
pling is exponentially decaying, rather than just nearest 
neighbour. The encouraging numerical results we obtain 
motivate us to present a possible extension of the quasi- 
diagonal approach to continuous space-time extended dy- 



namical systems (i.e. partial differential equations). 

For the convenience of the reader, we present our ap- 
proach in a largely self-contained manner, leading in a 
natural progression from a scalar method, which gives 
very poor results, to a successful quasi-diagonal estima- 
tion of the Jacobian from multivariate time series. In the 
process we highlight the relationship between the "global- 
ly" of the fitting basis and the "curse of dimensionality" . 
We also emphasise the importance of properly addressing 
the effects of the boundaries of the subsystem where we 
collect the multivariate time series. 

The paper is organised as follows. The next section 
gives a general introduction to spatio-temporal systems 
and describes the prototype model (a coupled map lat- 
tice) that is used in our numerical investigations. In sec- 
tion III we provide a short overview of Lyapunov spectra 



for extended dynamical systems and their relationship 
to fractal dimensions and the Kolmogorov-Sinai entropy. 
We give some examples and explain how to estimate the 
Lyapunov spectrum from subsystem information using a 
suitable rescaling. In the following section (IV) we at- 
tempt to estimate the Lyapunov spectrum using a scalar 
time series. We find that the lack of spatial information 
hinders any attempt to obtain a meaningful reconstruc- 
tion of the dynamics, and hence to estimate the spec- 
trum. In section ^ we present a systematic numerical 
study of estimates of the Lyapunov spectrum using dif- 
ferent spatio-temporal reconstructions. We discover the 
importance of including spatial information in order to 
obtain reasonable reconstructions of the dynamics. Wc 
also point out that boundary effects on the measured sub- 
system have to be properly addressed by truncating the 
outer layer(s) of the estimated Jacobian. In section VI 
we turn our attention to the effect of passing from a local 
fit to a global fit when increasing the embedding dimen- 
sion of the spatio-temporal reconstruction. We conclude 
that the "curse of dimensionality" precludes any hope of 
a usable reconstruction when our fitting basis does not 
give a good global a pprox imation to the dynamics. We 
then go on in section VII to exploit the sparse structure 
of the Jacobian when the coupling is short range to de- 
velop a quasi-diagonal estimation technique for the Jaco- 
bian. Additionally, we use the spatial homogeneity of the 
system to dramatically increase the amount of effective 
data points available when performing a local fit. This 
circumvents both the "curse of dimensionality" and the 
error induced by not using an appropriate global basis. 
Section VIII is devoted to the generalisation of our ap- 
proach to systems that have non-local, but exponentially 
decaying coupling. Finally, in the last section, we pro- 
pose a natural extension of our method to the estimation 
of the Lyapunov spectra of partial differential equations. 
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II. SPATIALLY EXTENDED SYSTEMS 

The occurrence of chaos in spatio-temporal systems 
has recently attracted the attention of a large part of 
the dynamical systems community. There exists nowa- 
days a broad understanding of low-dimensional chaotic 
systems. However, the same cannot be said of high- 
dimensional systems and in particular of spatially ex- 
tended systems. The addition of a spatial extent to 
the dynamics produces a complex interplay between the 
local dynamics (the original dynamics before including 
spatial interactions) and the spatial interactions. Some- 
times, this interplay triggers the so called phenomenon 
of spatio-temporal chaos. Loosely speaking, this refers to 
systems that combine a familiar temporal chaotic evolu- 
tion with an additional decay of spatial correlations. In 
a spatio-temporal chaotic regime both space translations 
and time evolutions exhibit instabilities and it is even 
possible to define spatial and temporal Lyapunov expo- 
nents ||. One possible simple mechanism for obtaining 
spatio-temporal chaotic motion is to spatially couple low- 
dimensional chaotic units; although one has to be care- 
ful because the coupling sometimes tends to reduce the 
spatio-temporal instabilities 0. However, it is also possi- 
ble to produce spatio-temporal chaos through the spatial 
interaction of well-behaved (non-chaotic) units. This is 
the case for some metapopulation dynamics models [§) . 

In contrast with non-spatially distributed systems, 
spatio-temporal systems possess a spatial extent. This 
may be discrete, giving a lattice with a local dynamical 
unit at each site, or continuous. The typical model in the 
latter case (if time is also continuous) is a partial differ- 
ential equation (PDE). By discretising space we obtain 
a lattice of ordinary differential equations (commonly re- 
ferred to as a lattice differential equation). In the discrete 
space case, the system can be viewed as a collection of 
low dimensional dynamical systems coupled together via 
some spatial rule (for a review of current research see 
Ref. ji| ) . Examples of this kind of model are widespread 
in the literature, particularly in the field of solid state 
physics where they are used to study the dynamics of 
interacting atoms arranged in a lattice (see Ref. (!(]] and 
references therein). 

In this paper we focus our attention on a third category 
of extended dynamical systems where not only space but 
also time is discrete. In such a case the model consists 
of low-dimensional dynamical units with discrete time 
(i.e. maps) arranged in some discrete lattice configura- 
tion in one or more spatial dimensions. Such models 
are usually called coupled map lattices (CMLs). They 
were first introduced in 1984 as simple models for spatio- 
temporal complexity [pl]-p^[. Despite their computa- 
tional simplicity, CMLs are able to reproduce a wide 
variety of spatio-temporal behaviour, such as intermit- 
tency Q , turbulence and pattern formation (jl7| , 



to name just a few. 

In this paper we shall focus on a one-dimensional array 
of sites. The length of this array could be infinite but here 
we restrict ourselves to an array of length N with periodic 
boundary conditions. At the i-th site we introduce a 
discrete time local dynamical system whose state at time 
n we denote by a;™. We suppose that the same local 
map / acts at every spatial location (so that the local 
dynamics is homogeneous) . In the simplest case, the local 
variable x™ is taken to be one-dimensional. The dynamics 
of the CML is then a combination of the local dynamics 
and the coupling, which consists of a weighted sum over 
some spatial neighbourhood. The time evolution of the 
i-th variable is thus given by 

where the range of summation defines the neighbour- 
hood. The coupling parameters Sk are site-independent, 
and satisfy ^2 £ k = 1- The commonest choice for the 
coupling scheme is 

x? +1 = (1 - e )/(s?) + | + /OrJVO) , (2) 

which is sometimes called a diffusive CML. This is a dis- 
crete analogue of a reaction-diffusion equation. There is 
now a single coupling parameter e which is constrained 
by the inequality < e < 1, to ensure that the signs of 
the coupling coefficients in (||) (i.e. e/2 and 1— e) remain 
positive. 

III. LYAPUNOV SPECTRA FOR EXTENDED 
DYNAMICAL SYSTEMS 

Let us now give an overview on the extraction and 
application of Lyapunov spectra for spatio-temporal sys- 
tems. For a iV-dimensional dynamical system there exist 
N Lyapunov exponents which in principle can be ob- 
tained from the eigenvalues of the matrix 

T= lim \P(nf ■ P(n)] 1/2n , (3) 

n — >oo 

where P(n) corresponds to the product of the first n Ja- 
cobians along the orbit and ( ■ ) tr denotes matrix trans- 
pose. It is well known that given an invariant measure 
(which is usually assumed to be the natural measure 
for the dynamics) the Multiplicative Ergodic Theorem 
(e.g. Ref. [lq] ) ensures that the limit exists for almost all 
initial conditions x°, and if the measure is ergodic then its 
eigenvalues are independent of x° . The Lyapunov expo- 
nents are then defined as the logarithms of these eigenval- 
ues (which are obviously non- negative) . Although com- 
puting the exponents using equation (^) appears straight- 
forward, in practice multiplying the Jacobians is an ill- 
conditioned procedure since the most expanding direc- 
tion swamps all the other expansion/contraction rates. 
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Therefore one has to resort to algorithms that regularly 
reorthogonalise the product of the Jacobians. Such al- 
gorithms perform a QR decomposition every few iterates 
and are closely related to standard methods of computing 
eigenvalues The QR decomposition can be carried 
out using Gram-Schmidt, Householder or Givens based 
techniques. In this paper we use an efficient Householder 
method where only the terms required for the computa- 
tion of the Lyapunov exponents are actually calculated 

We define the Lyapunov spectrum (LS) as the set of 
Lyapunov exponents {A,}^ arranged in decreasing or- 
der. The LS not only gives the expansion/contraction 
rates of infinitesimal perturbations, but can also provide 
estimates of fractal dimensions and entropies. Thus, for 
instance, the dimension of the chaotic attractor (i.e. in- 
formally the effective number of degrees of freedom) is 
given by the Kaplan- Yorke conjecture through the 
Lyapunov dimension 



Dl = j 



1 



5> 



(4) 



where j is the largest integer for which Y^i=i ^ - > 0- 
It is also possible to extract an upper bound for the 
Kolmogorov-Sinai (KS) entropy h from the LS by the 
following approximation jlq] 



E 



a; 



(5) 



where the summation is over the positive Lyapunov ex- 
ponents A^~ . The KS entropy quantifies the mean rate of 
information production in a system, or alternatively the 
mean rate of growth of uncertainty due to infinitesimal 
perturbations. 

As an example, figure ^.a shows the LS for a coupled 
logistic lattice of size TV = 20 with periodic boundary 
conditions obtained using the known dynamics (^). In 
figure ^.b we plot the sum of the Lyapunov exponents 
y^l— i Afe from where the KS entropy can be estimated 
(maximum of the curve) and the Lyapunov dimension 
can be extracted (intersection with the horizontal axis). 
Notice that the Lyapunov dimension ~ 15.4 is com- 
parable to the total size of the lattice N = 20. Also 
observe that Dl ~ 15.4 implies the presence of a high 
dimensional attractor. This point will be addressed in 
more detail in the following sections. 

The computation of the whole LS for spatio-temporal 
systems is a cumbersome task, particularly as the size 
of the system increases. Even the most efficient meth- 
ods involve 0(N 3 ) arithmetic operations per time step 
Jl9| for a lattice of size N. Even for moderate lattice 
sizes (N ~ 50) the number of operations required is al- 
ready considerable. Moreover, taking into account that 
sometimes the (time) convergence of the Lyapunov ex- 
ponents is extremely slow, the time necessary to com- 
pute the whole LS for spatio-temporal systems quickly 



becomes prohibitive. Furthermore, one is often inter- 
ested in the behaviour of the system in the thermody- 
namic limit where the number of lattice sites N goes to 
infinity. In such a case it quickly becomes impossible to 
compute the whole LS and one has to resort to subsystem 
rescaling techniques that we shall now briefly describe. 
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FIG. 1. a) Lyapunov spectrum for a CML of N = 20 fully 
chaotic logistic maps f(x) — 4x(l — x) with coupling strength 
e — 0.4. b) The Lyapunov dimension Dl and the KS en- 
tropy h may be estimated from the sum of the Lyapunov 
exponents. In this case the largest Lyapunov exponent is 
Ai ~ 0.36, Dl ^ 15.4 and h ~ 1.8. 

Consider a one-dimensional lattice of coupled one- 
dimensional dynamical units. High dimensional cases 
may be treated in a similar way. Assume the array is 
large (N ^> 1) and suppose initially that it is made up of 
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smaller independent (uncoupled) subsystems of size N s . 
For simplicity take N to be a multiple of N s . The LS of 
the whole lattice then consists of N/N s exact copies (as- 
suming spatial homogeneity) of the subsystem LS of N s 
Lyapunov exponents. Thus, the LS for the whole lattice 
consists of N s Lyapunov exponents, each with multiplic- 
ity N/Ng. As we relax the assumption that the subsys- 
tems are uncoupled we hope that the LS does not change 
significantly p2p5| ]. If this is the case, the whole LS of 
the original lattice may be approximated by a rescaled 
version of its subsystem LS. 

An obvious choice for such a rescaling is to multiply the 
index i of the iV s -subsystem Lyapunov exponents Xi(N s ) 
by r' = N/N s . However, careful analysis of the LS of 
a homogeneous state suggests that a better rescaling for 
one-dimensional lattices is given by [M 



v n = (<p n ,<p n - 



-(<*- 



(7) 



N + l 



1 



(6) 



The rescaling of the whole LS from subsystem infor- 
mation leads us naturally to define intensive quantities 
from extensive quantities. As a direct consequence of 
this rescaling, using the same argument as in the previous 
paragraph, it is straightforward to see that the Lyapunov 
dimension and the KS entropy are extensive quantities 
(i.e. Dl and h increase linearly with N). Moreover, it 
is useful to introduce their respective densities by sim- 
ply dividing them by the system volume (lattice size) 
p5j . This even allows such quantities to be defined in 
the thermodynamic limit, although care has to be taken 
in doing this for systems such as PDEs where space is 
continuous p6| . 

By using subsystem rescaling one can considerably re- 
duce the computing resources required to estimate the 
whole LS of a large spatio-temporal system. Further- 
more, by restricting oneself to observing data from a sub- 
system, one reduces the dimensionality of the resulting 
time series and increases the likelihood of success in ap- 
plying reconstruction techniques. 



IV. LYAPUNOV SPECTRA FROM UNIVARIATE 
TIME SERIES 

In the previous section we assumed that the dynamics 
of the system whose LS we wished to compute was known. 
This is not the case for many complex physical systems. 
We therefore now turn to the task of estimating the LS 
when the only information available about the system is 
a time series of observed data. Let us start by apply- 
ing standard delay reconstruction techniques for a single 
observable time-series. Thus suppose that the univariate 
time series {</?"} has been measured from a k dimensional 
system. We can create a high dimensional reconstructed 
state space by constructing the delay vectors 



where d is the so-called embedding dimension. Takens' 
Theorem J27J says that for d > 2k + 1, the dynamics 
induced on such delay vectors is generically smoothly 
conjugate to the original dynamics. This apparently ar- 
bitrary condition is a simple geometrical disentangling 
in order to avoid self intersections on the reconstructed 
manifold p8| . The smooth conjugacy between induced 
and original dynamics asserts that the time series con- 
tains all the coordinate-free properties of a dynamical 
system. Specifically, since the LS is invariant under 
smooth conjugacy (i.e. smooth coordinate changes) the 
LS computed from the dynamics of the delay vectors will 
be the same as that of the original system. Note that the 
time series need not to be the time history of one of the 
system's state variables but can be any generic function 
of the state of the system. In particular, in principle it is 
possible to reconstruct all the Lyapunov exponents from 
a univariate time series. To anyone familiar with the 
proof of Takens' Theorem is it obvious that it generalises 
straightforwardly to multi-variate time series, though as 
far as we are aware, a proof has never appeared in the 
literature. 

Unfortunately, Takens' Theorem does not actually ap- 
ply to the CMLs considered here Q. Firstly, the maps 
that we use as local dynamics are not invertible, and 
the proof of Takens' Theorem depends fundamentally on 
such invertibility. Secondly, we assume translation in- 
variance of the CML (i.e. the dynamics and coupling are 
spatially homogeneous) and such symmetry renders our 
systems far from generic. Finally, in a similar fashion, 
measurement functions depending on only a single site 
(or group of sites) are not generic in the space of all ob- 
servables, indeed a generic function will depend on all N 
sites. On the other hand our results, and in particular 
the comparison with the LS computed in the previous 
section assuming full knowledge of the actual dynamics, 
suggest that we are able to estimate the real LS from 
observed time series. We shall therefore use delay recon- 
struction techniques throughout this paper without wor- 
rying about this failure of Takens' Theorem (and in any 
case the theorem only draws conclusions for generic dy- 
namics and measurement functions, and so cannot guar- 
antee reconstruction for any particular time series). 

For the extended dynamical systems we are consider- 
ing, the dimension k is the size of the lattice N times 
the dimension of the local dynamics. Hence for the lo- 
gistic coupled map lattice in figure [l] we get k = N. 
Even for quite moderate N an embedding dimension of 
2N + 1 is clearly quite impractical. Of course, the con- 
dition d > 2k + 1 is only a sufficient one, and it is well 
known that some systems (e.g. the Lorenz equations) can 
be reconstructed for smaller values of d. Furthermore, 
Sauer and Yorke have a sharper estimate for d that 
is still sufficient to preserve the box-counting dimension 
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and the Lyapunov exponents 

d>D B +D T , (8) 

where Db is the box-counting dimension of the attrac- 
tor and Dt is a "tangent dimension" which is informally 
the maximum (over all points) dimension of the tan- 
gent space of the underlying dynamics. However even 
this is still very prohibitive. For example, for the sys- 
tem in figure |lj) the Lyapunov dimension for even quite 
a small lattice of N = 20 coupled logistic maps we have 
Dl ~ 15.4. Then assuming Db = L>l [i-e. essentially the 
original Kaplan- Yorke conjecture), the bound (|J) gives 
d > 15.4 + Dt, and one would typically expect Dt to 
lie between Db and N, giving a value of d from 30 up- 
wards. Of course the Sauer- Yorke inequality is also just a 
sufficient condition, but if one is to preserve all N expo- 
nents, it seems that an absolute minimum of d > N will 
be necessary, and even this is too large to be generally 
practicable. 
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FIG. 2. Estimates of the LS for a logistic coupled map 
lattice of size N = 20 and coupling e = 0.4 computed using a 
time delay reconstruction (with d € [32, 42]) from a univariate 
time series (thin lines) . A second order local fit from a sample 
of 10 4 points was used and the number of neighbours was 
set to the number of parameters in the fitting process plus 
20. For comparison the spectrum computed from the original 
dynamics is depicted with a thick line. An orbit length of 10 4 
iterations was used in both cases and the computed LS was 
rescaled using (^). 

Indeed, as indicated in the introduction, an embed- 
ding dimension much larger than d ~ 6 will lead to the 
so-called "curse of dimensionality" . Thus in order to esti- 
mate the LS from the time series we need to approximate 
the map that governs the dynamics of the delay vectors, 
that is the map F such that F(y n ) = y n+1 . There ex- 
ist several methods for doing this, of which the best for 



the purpose at hand are all based on local approxima- 
tions. These can be constructed by scanning the time 
series to find a set of past close encounters (correspond- 
ing to neighbours in the reconstructed state space) of 
the current reconstructed state y n . Once these neigh- 
bours have been found, an approximation F to F at y n 
can be computed using, for example, a least squares fit 
within an appropriate space of functions (e.g. linear or 
quadratic). For a comprehensive review see the recent 
book by Kantz and Schreiber J3(|. For large embedding 
dimensions (d ~ 6) it becomes extremely difficult to find 
close neighbours and hence, due to data limitations, the 
approximation to F ceases to be local. This in turn can 
lead to large errors in the determination of the Lyapunov 
exponents. 

It is clear then that function approximation techniques 
are bound to fail when reconstructing spatio-temporal 
systems where the dynamics is typically high dimen- 
sional. As an example figure |^ shows the results of an 
attempt to estimate the whole LS from a univariate time 
series obtained from a single site of a couple logistic lat- 
tice. The measurement function in this case is the pro- 
jection onto a single site, so that <p n = x™ for a fixed j. 
An embedding dimension in the range 32 < d < 42 was 
used. It is obvious from the figure that this temporal 
delay reconstruction of the LS (overlapping thin lines) 
completely fails to reproduce the LS of the original sys- 
tem (thick line). Other values of the reconstruction pa- 
rameters than the ones given in the figure caption were 
also tried yielding the same qualitative results. 

V. LYAPUNOV SPECTRA FROM 
SPATIO-TEMPORAL DATA 

Since the direct application of univariate (temporal) 
delay reconstructions clearly fails for spatio-temporal sys- 
tems, let us now try to exploit the spatial extent of the 
system and use a spatio-temporal delay reconstruction. 
The framework for this is presented in Q . Here, we shall 
show the benefits of including spatial information in esti- 
mating the LS. We shall also investigate the advantages 
of truncating the outer layer of a reconstructed Jacobian 
in order to avoid boundary effects. We have already ad- 
dressed these issues in a previous paper [Q], however we 
repeat some of the analysis here for completeness. 

A spatio-temporal reconstruction is obtained by re- 
placing the delay vectors (0) of the previous section by 
the spatio-temporal delay vectors 

V? = (#,#_!,..., #_ (< «._i)), (9) 

whose entries 0" = (x™, a;™ -1 , . . . , x™ ^ dt ) are time- 
delay vectors and the spatial index i is fixed. The overall 
embedding dimension for such a spatio-temporal recon- 
struction is d = d s dt, where d s and dt denote the spatial 
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and temporal embedding dimensions respectively. The 
standard temporal delay reconstruction (Q) is recovered 
by setting d a = 1. Spatio-temporal reconstructions have 
proved useful in various contexts ]3l|-|33|] . By adding the 
extra spatial delay in the time series one is including im- 
portant information that is otherwise very difficult to ob- 
tain. In fact, it is extremely difficult to extract informa- 
tion about neighbouring dynamics just by doing spatially 
localised measurements. This is because of the rapid de- 
cay of spatial correlations in locally coupled extended 
dynamical systems. In a previous paper || we showed 
how this extremely rapid decay in correlations allows a 
quite small truncated lattice with random inputs at the 
boundaries to reproduce the local dynamics of a large, 
potentially infinite, lattice. This tends to suggest that a 
univariate time series cannot in practice reconstruct the 
dynamics of a whole spatio-temporal system. 
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FIG. 3. LS estimated from spatio-temporal delay recon- 
structions for a logistic coupled map lattice. The LS com- 
puted from the original dynamics is depicted with a thick line 
while the different spatio-temporal delay reconstructions are 
depicted with dashed lines. All the fitting parameters are 
the same as in figure ^| (second order fitting, 10 4 points, 10 4 
iterates) . 

In figure || we present several estimates of the LS of 
the coupled logistic lattice (as in figure |J) for different 
spatio-temporal embedding dimensions d. The different 
reconstructions correspond to a broad sample of choices 
of d s and d t . Note that simply increasing d = d s dt does 
not improve the estimates. In fact, in order to get a good 
estimate of the LS, it seems preferable to take dt = 1 and 
d s as large as possible. This is indicated by the empty 
and filled circles corresponding to dt = 1 and d s = 10 
and d s = 20 respectively. In the figure we show only a 
few combinations of (d s , dt). However, we originally con- 



sidered a much broader choice without finding any qual- 
itative differences. We therefore conclude that a spatio- 
temporal reconstruction with dt = 1 (that is a pure spa- 
tial delay reconstruction) gives the best results. This is 
easy to explain since we are using the actual dynamical 
variables of the system (ip™ — a;™) as observables and 
are measuring them in a complete window. Thus, we do 
not need to increase the embedding dimension to avoid 
self intersection in the reconstructed space since we are 
already in a natural space for the system. 
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FIG. 4. LS estimated from pure spatial delay reconstruc- 
tions for a logistic coupled map lattice. All the fitting pa- 
rameters are the same as in figure || (second order fitting, 10 4 
points, 10 4 iterates). 

In figure ^ we show the reconstruction of the LS using 
a pure spatial reconstruction with increasing d s . These 
estimates are much more promising than before. Note 
that since our original system is of size N = 20 it is 
not possible to choose d s > 20. Indeed, we have to be 
careful, since taking d t — 1 and d s = N = 20 corresponds 
to measuring the whole state space and is therefore not 
a genuine delay reconstruction. This explains why the 
reconstructed LS for dt = 1 and d s — N — 20 agrees so 
well with the direct computations (triangles in figure ^) 
though note that in figure [l] we know the Jacobian of the 
dynamics exactly, whilst here we still have to do some 
function fitting. 

From figure || we see that as d s is increased a better 
estimate of the LS is obtained. One may think that this 
is simply the effect of getting a higher embedding dimen- 
sion. However, the reason for the apparent convergence 
of the reconstructed LS towards the exact LS as d s — > N 
is simply due to a reduction of the noise from boundary 
effects. Consider the sites at the boundary of our spatial 
window in which we measure the multivariate time se- 
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ries. The dynamics of these sites depends on sites whose 
evolution is not explicitly contained in the measurement 
window. As a consequence, the estimates of the Jacobian 
for these sites contain spurious terms and are subject to 
error. For instance, the fitting procedure for the first row 
tries to compensate for the lack of information from the 
left neighbouring site that is outside of the measurement 
window by including an artificial dependence on all sites 
inside the window. This will yield erroneous non-zero 
terms off the tri-diagonal. 
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FIG. 5. LS estimated from pure spatial delay reconstruc- 
tions for a logistic coupled map lattice using truncation of 
the outer layer of the Jacobian. All the fitting parameters are 



the same as in figure 
iterates) . 



(second order fitting, 10 points, 10 



We can avoid these effects by observing that the 
boundary sites contribute towards the Jacobian at its 
outer rows and columns. We can therefore remove such 
boundary effects by simply truncating the outer rows 
and columns of the reconstructed Jacobian to yield a 
(d s — 2) x (d s — 2) matrix Note that one must first 
fit the d s x d s Jacobian and only then truncate its outer 
rows and columns. Also observe that for a spatial delay 
embedding dimension d s we obtain d s —2 Lyapunov expo- 
nents. The results of applying this technique arc shown 
in figure |^. Comparing this to figure |] where the outer 
layer of the Jacobian was not truncated we see a dra- 
matic improvement, and in particular can obtain good 
estimates even for quite small d s (see also Ref. &), It 
is obvious that the truncation of the outer layer of the 
Jacobian has to coincide with the range of the coupling: 
the larger the coupling range the more outer rows and 
columns must be truncated. This is considered further 



below in section VIII 



VI. LOCAL VERSUS GLOBAL FITTING: A 
MATTER OF DIMENSIONALITY 

It is somewhat surprising that the estimated LS in fig- 
ure H are in such good agreement with the actual LS 
despite the fact that the embedding dimensions are a) 
too small to disentangle the state space and b) too large 
to avoid the "curse of dimensionality" . The explanation 
for this apparent contradiction is that the measurement 
functions are the actual variables of the system and the 
form of the dynamical equations is actually contained in 
the space spanned by the basis functions used to approx- 
imate them. In particular, the dynamics of our CML is 
quadratic in nature (c.f. equation (^) with logistic local 
maps) and we are using a quadratic fit. Thus, a local fit 
is in fact also a global one and the problem of not finding 
enough neighbours in high-dimensions is avoided. It thus 
turns out to be possible to obtain very good estimates of 
the dynamics with even moderate amounts of data. 
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FIG. 6. Skewed logistic map (thick line) for for b = 0.2 
obtained by applying ( |Io| ) to the fully chaotic logistic map 
(dashed line). 

This situation is however unlikely to arise in most real 
applications, and typically we expect that the (delay re- 
constructed) dynamics will not lie in the space spanned 
by the basis functions. To investigate the effects of this 
we now turn our attention to a skewed logistic local map 
which yields a CML whose dynamics is not contained in 
the space of quadratic polynomial functions that we use 
for fitting. This map is obtained from the standard logis- 
tic map f{x) = ax{\ — x) by applying the following skew 
transformation of the unit square 



K(x,y) = (x + by,y), 



(10) 



where & is a parameter determining the degree of skew. 
In this way we obtain the map (figure ||) 
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9b (x) 



-1 + ab(2x - 1) + + ab) 2 - 4abx 
2a~& ' 



(11) 



where b must satisfy — 1/a < b < 1/a so that the deriva- 
tives of gb at and 1 remain bounded. 
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FIG. 7. Estimation of the LS for a coupled map lattice 
of skewed logistic maps. The size of the lattice is N = 20, 
the coupling parameter is e — 0.4 and the skew parameter is 
b = 0.2. The original LS is indicated by the solid line and the 
estimates using a range of spatial embedding dimensions d a 
are depicted by various symbols. For large d s the estimates 
deteriorate because of the "curse of dimensionality" . All the 
fitting parameters are the same as in figure ^ (second order 
fitting, 10 4 points, 10 4 iterates). 

We now repeat the estimation of the LS for this skewed 
map. Figure shows the effects of increasing the spatial 
embedding dimension d s in the case of large skew b = 0.2. 
Observe that for small d s (< 6) the approximation of the 
LS is quite good (see square points in figure). However, 
as d s increases, the estimate deteriorates. This reflects 
the fact that as the embedding dimension grows the fit- 
ting algorithm has more and more difficulty finding close 
neighbours. Since the skewed logistic map cannot be well 
approximated globally by a quadratic function the LS es- 
timate deteriorates. Note however that even though an 
accurate estimate of the LS is not possible in this case, 
the estimated LS does have the right qualitative features. 
This is due to the fact that the skewed logistic map still 
resembles a quadratic function to some extent and thus 
a global fit is still reasonable. This is illustrated in fig- 
ure H where we present the effect of varying the degree 
of skew 6, and hence the degree of departure from the 
space of quadratic functions. We see that as the skew 
is increased, the error between the actual and the esti- 
mated LS rapidly increases. Using a local map which is 
even further from the space of quadratic ones leads to yet 



greater discrepancies [pi. 

Finally recall also that the number of Lyapunov expo- 
nents obtained for a given d s is d s — 2 due to the trunca- 
tion of the outer layers of the Jacobian. Therefore if we 
try an avoid the "curse of dimensionality" by working in 
low embedding dimensions with d s < 6 we obtain only 
4 points or less on the LS density curve. Even if these 
are accurate, we cannot in general expect to extrapolate 
the whole spectrum from just these 4 points. Hence esti- 
mates of Lyapunov dimension and KS entropy are likely 
to be quite poor 0. 
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FIG. 8. Effect of the skew of the local map when estimat- 
ing the LS. Curves correspond to the LS computed with the 
known dynamics and the symbols depict various pure spatial 
delay reconstructions. All the fitting parameters are the same 
as in figure^ (second order fitting, 10 4 points, 10 4 iterates). 

The results presented in this section suggest that as 
soon as the dynamics cannot be globally approximated 
by the space spanned by the basis functions used in the 
fit, and as soon as more than a handful of Lyapunov 
exponents is required, standard spatio-temporal embed- 
ding techniques cannot reliably estimate the LS. In the 
next section we shall show how this difficulty can be over- 
come by focusing on the estimation of only the nontrivial 
entries in the Jacobian. An alternative approach is to 
try and find an adapted fitting basis that will allow us 
to fit the dynamics globally. This can be done by first 
extracting an estimate of the local dynamics from the 
time series and then using this to construct an appro- 
priate set of basis functions. One possible way of esti- 
mating the local dynamics is to use time-delay plots [||. 
An even more promising technique is to consider quasi- 
homogeneous states in a small window and their time 
evolution 0. 
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VII. QUASI-DIAGONAL RECONSTRUCTION OF 
THE JACOBIAN 



In this section we improve the method presented above 
by making use of the local nature of the coupling in our 
CML, which results in the Jacobian exhibiting a banded- 
diagonal structure. It is thus natural to only attempt to 
estimate the non-zero entries in this Jacobian. We call 
this a quasi- diagonal reconstruction. This allows us to 
carry out the local fit in a low dimensional space, and 
avoids the difficulties described above. A similar idea 
was employed by Biinner and Hegger H. However, they 
only applied it to a standard logistic CML, which as we 
have seen above is a relatively easy system for which to 
estimate the spectrum using conventional means. It is 
thus impossible to judge from their work how much of a 
benefit the quasi-diagonal approach actually brings. Ad- 
ditionally, here we combine this technique with a trun- 
cation of the outer layers of the Jacobian, which, as we 
have seen above, brings substantial benefits, but was not 
considered in 

The idea behind a quasi-diagonal reconstruction of the 
Jacobian is quite simple. The entries of the Jacobian 
J(n) at time n are given by 



Jki{n) 



(12) 



Due to our assumption of localised coupling most of these 
terms are identically zero. The only non-zero terms lie 
within a distance of the diagonal given by the distance 
over which the coupling acts. In the particular case of a 
nearest neighbour CML (||), the Jacobian is tri-diagonal, 
and has only three non-zero elements in each row: 



Jin) 
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dt 1 



V o 



d 
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d., 



(13) 



where, for simplicity, we use the notation: 

QU 



8r n+1 



ux i +v-l 

Note that the Jacobian (|l3| ) is extracted from sites io to 
io + d s — 1 and so is of size d s x d s . The starting index ig 
is arbitrary since we are assuming spatial homogeneity. 

We estimate this Jacobian in a row by row fashion by 
fitting local dynamics of the form 



for (1 < r < d s ). This fit is carried out in just a three 
dimensional space, as opposed to a d s dimensional as be- 
fore, and close neighbours are defined by their distance 
from (x^-i) x r> x r+i)- This thus avoids any problems 
with high-dimensionality. Furthermore, if the dynamics 
is translation invariant (i.e. spatially homogeneous), we 



can use any triple of the form (x r : 



i+i), regard- 



less of position. This dramatically increases the effective 
amount of data available to us for performing the fit. Ad- 
ditionally, if the dynamics is also isotropic, i.e. invariant 
under the symmetry that exchanges x" with x r ^ [ _-, we 
can use all triples in reverse order (a^ +1 , x 7 -, Of 
course, if the coupling acts over a longer range then the 
map F r will depend on more variables, but as long as the 
coupling remains reasonably local this approach will still 
bring advantages (we investigate this further in the next 
section). 



-0.2 




(14) 
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FIG. 9. Effect of truncating the outer row of the Jacobian 
for a tri-diagonal estimate of the LS for a lattice of coupled 
skewed logistic maps as in figures [?] and ^. The skew is b = 0.2 
and the width of the window used to extract the multivariate 
time series is d„ sites. The open symbols correspond to using 
the whole fitted Jacobian while the solid symbols correspond 
to truncating the outer layer of the Jacobian. All the fitting 
parameters are the same as in figure ^| (second order fitting, 
10 4 points, 10 4 iterates). 

The first and last rows of J(n) contain only two non- 
zero entries due to the lack of information from outside 
the observation window. As before, the entries will be 
estimated incorrectly since the fitting algorithm will try 
and compensate for the lack of the missing entry by erro- 
neously adjusting the remaining two. We therefore trun- 
cate the outer layers of the Jacobian as in the previous 
section and compute only d s — 2 Lyapunov exponents 
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from a window of size d s . The improvements due to such 
truncation are shown in figure ^. Note that these be- 
come less apparent as the width d s of the observation 
window increases. This is consistent with the fact that 
the boundary effects become less significant as the size 
of the subsystem grows M. Nevertheless, given that it 
always leads to better estimates, we shall continue to 
employ such truncation throughout the remainder of the 
paper. 
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FIG. 10. Tri-diagonal estimate of the LS for a lattice of 
coupled skewed logistic maps as in figures [j] and || The skew 
is b — 0.2 and the width of the window used to extract the 
multivariate time series is 20 sites. All the fitting parameters 
are the same as in figure ^| (second order fitting, 10 4 points). 

Figure |l^ shows an example of the estimation of the 
LS using a quasi-diagonal reconstruction with a large 
window d s = 20. Compared to figures ^ and |8| we see 
that this approach can give an excellent estimate, even 
though the local map cannot be well approximated by a 
quadratic fit. Similar results where obtained for different 
maps and coupling strengths. 



VIII. ESTIMATING THE LYAPUNOV 
SPECTRUM FOR EXPONENTIALLY DECAYING 
COUPLING 

The CML used to illustrate our method in the pre- 
vious sections only allowed interactions between nearest 
neighbours. There are however many other systems of 
interest where the coupling acts over greater distances. 
In this section we investigate the performance of our ap- 
proach in such cases. Typically, it is assumed that the 
coupling decreases in strength with distance, giving a 
band-diagonal Jacobian with sub-diagonals whose entries 



decay as we move away from the diagonal. A suitable 
paradigm model to represent this is a CML where the 
coupling range is in fact infinite but whose coupling co- 
efficients decay exponentially with distance: 



1-/3 
1 + 



£ |i_fc| /W-*). 



(15) 



where the coupling parameter /3 £ (0, 1). The limit (3 — > 
corresponds to the uncoupled case, i.e. depends 
only on x™_ k . The limit j3 — > 1 corresponds to global 
coupling of all the sites with the same coefficient. Thus 
increasing [3 effectively increases the range of the cou- 
pling and results in more and more sub-diagonals of the 
Jacobian becoming significant. Whilst for (3 close to 0, 
it is reasonable to expect that a tri-diagonal reconstruc- 
tion will still suffice, it is clear that as (3 increases, more 
sub-diagonals will need to be taken into account. 

We thus investigated the effects of using a higher- 
diagonal reconstruction for data from a lattice with a 
relatively large (3 (0.35). If we fix the size q of the ad- 
missible interaction range we need to estimate a map 
depending on 2q + 1 variables: 



r K^r — q ' ' ' ) "^r— 1' r* ' r+1 ' ' ' ■ ' r+qv 



(16) 
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FIG. 11. Estimation of the LS using a (2q + l)-diagonal 
reconstruction of the Jacobian for a lattice of coupled skewed 
logistic maps with decaying exponential coupling (|l5|). The 
skew is b = 0.2 and the decay rate is (3 — 0.35. The cou- 
pling parameters for (3 = 0.35 correspond to eo — 0.481, 
£±i ~ 0.168, £± 2 ~ 0.059, £± 3 ~ 0.021, £ ±4 ~ 0.007, 
in the general CML formulation (|l]). All the fitting param- 
eters are the same as in figure [| (second order fitting, 10 4 
points, 10 4 points). 

We now need to discard q outer layers of the estimated 
Jacobian to remove boundary effects, and hence if our 
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observation window is of size d s , we can estimate d s — 2q 
Lyapunov exponents. 

Figure [ll] shows the results of this approach applied 
to the CML coupling scheme §U^) for different values of 
q (still with a skewed logistic map for local dynamics). 
We see that the tri-diagonal (q = 1) reconstruction com- 
pletely fails to approximate the LS. However, as q in- 
creases the reconstruction rapidly improves. For q = 2 
the largest Lyapunov exponents are captured accurately 
and the estimate of the remainder of the spectrum is still 
quite poor. This improves significantly for q = 3 and 
g = 4 where the approximation is rather good, particu- 
larly given that the coupling is not genuinely over a finite 
range. Increasing q further led to a deterioration of the 
estimates of the LS. This is due to the reappearance of 
the "curse of dimensionality" : for q > 5 we have to find 
neighbours in dimension <i=2g+l>ll. We repeated 
these experiments for other values of (3 and other local 
maps, obtaining similar qualitative results. 

IX. DISCUSSION AND GENERALISATIONS 

We have shown that by using a quasi-diagonal recon- 
struction of the Jacobian we are able to estimate the LS 
of a variety of CML's using only time series data with 
a reasonable degree of accuracy. The most difficult test 
was where the local dynamics was given by a skewed lo- 
gistic map and the coupling was exponentially decaying. 
In that case, we had to take a 9-diagonal reconstruction 
of the Jacobian (i.e. q = 4). Increasing q beyond this 
leads to increasing errors, due to the increase in the di- 
mensionality of the reconstruction space that we use. We 
therefore have a dichotomy when trying to estimate LS 
for systems with extended interactions: on the one hand 
we would like to include as many sub-diagonals as pos- 
sible (i.e. q as large as possible), however, on the other 
hand, q cannot be chosen too large (q < 5) otherwise we 
again encounter the "curse of dimensionality" . 

Of course, in practice if we do not know the dynamics 
of a system, we are unlikely to know a priori the range 
of the coupling. However, this is something that can 
be estimated from a multivariate time series, using for 
example some kind of cross correlation. This in turn will 
allow us to estimate the width of the band of the Jacobian 
that is most appropriate to capture the essence of these 
interactions. 

A related problem arises when trying to reconstruct 
LS from multivariate time series produced by systems 
with a continuous space variable, or when probing real 
life extended dynamical systems. Consider for example 
a one-dimensional PDE and assume that we are free to 
choose the sampling intervals in both time and space at 
which data is observed. The Jacobian will contain a dif- 
ferent number of significant sub-diagonals depending on 



the choice of these intervals. This is illustrated in fig- 
ure [l^ where we show the so called cone-horizon for the 
evolution of a perturbation in space-time. The cone- 
horizon corresponds to the region of space-time where 
a perturbation, applied at the summit of the cone, can 
influence downstream positions. Any point outside the 
cone-horizon cannot "feel" the presence of the distur- 
bance. Such cone-horizons always exist for extended dy- 
namical systems where the interactions have finite range, 
i.e. where information propagation has finite speed. Now 
suppose that we observe the system at regular intervals 
in space and time. This corresponds to the the intersec- 
tions of the grid shown in the figure. Let us denote by r 
the temporal sampling; for simplicity we assume the spa- 
tial sampling interval is fixed, since varying r is sufficient 
to demonstrate our point. For the example depicted in 
the figure a sampling interval of r yields a tri-diagonal 
Jacobian. This is because there are only three sites in 
the cone-horizon after a time r. Thus, for this choice 
of space-time discretisation it should be sufficient to use 
a tri-diagonal reconstruction of the Jacobian, i.e. with 
q = 1. However, if we double the sampling interval to 2r 
we shall need to use a 5-diagonal reconstruction in order 
to include the 5 sites in the cone-horizon after a time 2t. 

space 
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FIG. 12. Cone-horizon for a disturbance in a one-dimen- 
sional, space-time continuous, extended dynamical system. A 
disturbance applied at the summit of the cone may alter the 
downstream dynamics inside the cone-horizon (shaded area). 
Applying a given space-time discretisation (dashed mesh) in- 
duces a particular neighbourhood structure: a node at time 
t + t depends on a fixed number of nodes at time t (3 in 
this case). The slope of the cone- horizon fixes the number of 
sub-diagonals required in the quasi-diagonal estimation of the 
Jacobians. 

We intend to investigate the application of the quasi- 
diagonal method to the estimation of a LS from a time 
series generated by a PDE in a future paper. Addition- 
ally, we have hitherto assumed that the observable is the 
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actual state variable at a site. This is unlikely to be the 
case in many practical applications. In such a case, we 
expect that the best approach would be a combination 
of the quasi-diagonal reconstruction method presented 
here with a local spatio-temporal reconstruction of the 
dynamics. 
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